c2.unb=function(data,ybari,s2i,mi,Mi,M,n,N,param){
  data=as.name(data); N=N; mi=mi; n=n
  M=M; Mi=Mi; fpc=1-n/N; ybari=ybari
  s2i=s2i; Mbar=M/N
  muhat=(1/Mbar)*(sum(Mi*ybari))/n
  s2b=sum((Mi*ybari-Mbar*muhat)^2)/(n-1)
  if(param=='mean'){
    est=muhat
    vhat=fpc*s2b*(1/(n*Mbar^2))+(1/(n*N*Mbar^2))*sum(Mi^2*(1-mi/Mi)*(s2i/mi))
    bound=2*sqrt(vhat)
    lower=est-bound; upper=est+bound
  } else if(param=='total'){
    est=M*muhat
    vhat=M^2*(fpc*s2b*(1/(n*Mbar^2))+(1/(n*N*Mbar^2))*sum(Mi^2*(1-mi/Mi)*(s2i/mi)))
    bound=2*sqrt(vhat)
    lower=est-bound; upper=est+bound
  } else if(param=='proportion'){
    phati=yi/mi; qhati=1-phati
    est=sum(Mi*phati)/M
    s2r=sum((Mi^2*(phati-est)^2))/(n-1)
    vhat=fpc*s2r*(1/(n*Mbar^2))+(1/(n*N*Mbar^2))*(sum(Mi^2*(1-mi/Mi)*(phati*qhati/(mi-1))))
    bound=2*sqrt(vhat)
    l=est-bound; u=est+bound
    lower=sprintf("%.2f%%",l*100); upper=sprintf("%.2f%%",u*100)
  } 
  cat("","\n","Results from 2-stage Cluster design: Data =",data,"\n","\n","N =",N,"n =",n,"FPC =",fpc,"\n",
      "Estimate of",param,"=",est,"\n","Vhat(", param,") =",vhat,"\n","Bound =",bound,"\n","Lower Bound =",
      lower,"Upper Bound =",upper,"\n","")
  results=list(est=est,data=data,n=n,N=N,fpc=fpc,vhat=vhat,bound=bound,lower=lower,upper=upper,param=param,
               s2i=s2i,s2b=s2b)
}

# to use the function with its call:
# c2.unb(data,ybari,s2i,mi,Mi,M,n,N,param)

# data: name of dataset, in quotes
# ybari: numeric vector of the ith cluster means
# s2i: numeric vector of the ith cluster variances
# mi: numeric vector of the ith cluster elements
# Mi: numeric vector of the population cluster elements
# M: single number of total population elements
# N: single number of population size
# param: c('mean','total','proportion')